function  [D_CS, D_Ext, D_Mis_marginal, D_Mis_infra, D_Mis, D_SW, m_theta, m_D_E, m_D_E_marginal, m_D_E_infra, sd_theta, sd_D_E, sd_D_E_marginal, sd_D_E_infra, corr_theta_D_E, corr_theta_D_E_marginal, corr_theta_D_E_infra, D_Mis_SumStats, D_Mis_SumStats_1, D_Mis_SumStats_2, D_Mis_SumStats_1_marginal, D_Mis_SumStats_2_marginal,D_Mis_SumStats_1_infra, D_Mis_SumStats_2_infra ]  = Total_Welfare_decomp_WTP_EEgap_standard_marginal(weight,param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,MT,price, price_denom, price_standard_denom, rebate_estar_denom, estar_denom, kwh, kwh_standard, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom)
                                                                                                                                                                                                                                                    
    param_FE_k=param_FE;
	param_FEdemo_k=param_FEdemo;

    d_ext_J = repmat(d_ext',J,1);

dim_k=length(weight)
%parfor k=1:dim_k
for k=1:dim_k

	param_k=param{k};
    theta_k(k)=param_k(3);

diff_kwh = kwh_standard-kwh;
marginal_id_tmp = find(diff_kwh<0);
marginal_id=zeros(J,MT);
marginal_id(marginal_id_tmp)=1;
infra_id=zeros(J,MT);
infra_id    = 1-marginal_id;
out_id      = 0*marginal_id;
marginal_id = marginal_id';
infra_id    = infra_id';
out_id      = out_id';

%Welfare evaluated at the baseline 
[CS_0_i{k}, CS_0_avg(k), E_0_i{k}, E_0_avg(k), E_0_marginal_i{k}, E_0_infra_i{k}, E_0_out_i{k},  W_0_i{k}, W_0_avg(k), Error_0_i{k}, Error_0_avg(k), ext_0_i{k}, ext_0_avg(k), gvt_0_i{k}, gvt_0_avg(k), SW_0_i{k}, SW_0_avg(k), P_0_decision{k}, P_0_decision_j{k}, P_0_marginal_decision{k}, P_0_infra_decision{k}, P_0_out_decision{k}, share_marginal_0{k}, share_infra_0{k}, share_out_0{k}] = ...
Welfare_decomp_k_WTP_EEgap_standard(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, marginal_id, infra_id, out_id);

%Welfare evaluated at the standard
%kwh_standard
%price_standard_denom
[CS_i{k}, CS_avg(k), E_i{k}, E_avg(k), E_marginal_i{k}, E_infra_i{k}, E_out_i{k},  W_i{k}, W_avg(k), Error_i{k}, Error_avg(k), ext_i{k}, ext_avg(k), gvt_i{k}, gvt_avg(k), SW_i{k}, SW_avg(k), P_decision{k}, P_decision_j{k}, P_marginal_decision{k}, P_infra_decision{k}, P_out_decision{k}, share_marginal{k}, share_infra{k}, share_out{k}] = ...
Welfare_decomp_k_WTP_EEgap_standard(param_k,param_FE_k,param_FEdemo_k,id_ij,J,p_elec,t_pigou,d_ext_J,elec,lifetime,r,MCPF,MT,price, price_standard_denom, rebate_estar_denom, estar_denom, kwh_standard, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, marginal_id, infra_id, out_id);

E_0_ik(:,k)            = E_0_i{k};
E_ik(:,k)              = E_i{k};
D_CS_ik(:,k)           = CS_i{k}-CS_0_i{k};
D_Ext_ik(:,k)          = d_ext.*(E_i{k}-E_0_i{k});
D_Mis_marginal_ik(:,k) = (theta_k(k)-1)*elec.*(share_marginal{k}.*E_marginal_i{k}-share_marginal_0{k}.*E_0_marginal_i{k});
D_Mis_infra_ik(:,k)    = (theta_k(k)-1)*elec.*(share_infra{k}.*E_infra_i{k}-share_infra_0{k}.*E_0_infra_i{k});
D_Mis_ik(:,k)          = (theta_k(k)-1)*elec.*(E_i{k}-E_0_i{k});
D_SW_ik(:,k)           = SW_i{k}-SW_0_i{k};

%D_CS_avgk(k)   = CS_avg(k)-CS_0_avg(k);
%D_Ext_avgk(k)  = mean(mean(d_ext))*(E_avg(k)-E_0_avg(k));
%D_Mis1_avgk(k) = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k));
%D_Mis2_avgk(k) = theta_k(k)*t_pigou*E_avg(k);
%D_Mis_avgk(k)  = (theta_k(k)-1)*mean(elec)*(E_avg(k)-E_0_avg(k)) + theta_k(k)*t_pigou*E_avg(k);
%D_SW_avgk(k)   = SW_avg(k)-SW_0_avg(k);

D_E_marginal_ik(:,k) = share_marginal{k}.*E_marginal_i{k}-share_marginal_0{k}.*E_0_marginal_i{k};
D_E_infra_ik(:,k)    = share_infra{k}.*E_infra_i{k}-share_infra_0{k}.*E_0_infra_i{k};

end

 D_CS           = sum(weight'.*sum(D_CS_ik,1)/MT);
 D_Ext          = sum(weight'.*sum(D_Ext_ik,1)/MT);
 D_Mis_marginal = sum(weight'.*sum(D_Mis_marginal_ik,1)/MT);
 D_Mis_infra    = sum(weight'.*sum(D_Mis_infra_ik,1)/MT);
 D_Mis          = sum(weight'.*sum(D_Mis_ik,1)/MT);
 D_SW           = sum(weight'.*sum(D_SW_ik,1)/MT);

 E_0_k          = sum(E_0_ik,1)/MT;
 E_k            = sum(E_ik,1)/MT;
 D_E_marginal_k = sum(D_E_marginal_ik,1)/MT;
 D_E_infra_k    = sum(D_E_infra_ik,1)/MT;

 %D_CS_avg  = sum(weight'.*D_CS_avgk);
 %D_Ext_avg = sum(weight'.*D_Ext_avgk);
 %D_Mis1_avg = sum(weight'.*D_Mis1_avgk);
 %D_Mis2_avg = sum(weight'.*D_Mis2_avgk);
 %D_Mis_avg = sum(weight'.*D_Mis_avgk);
 % D_SW_avg  = sum(weight'.*D_SW_avgk);

%Compute sufficient stats
m_theta        = sum(weight'.*theta_k); 
m_D_E          = sum(weight'.*(E_k-E_0_k));
% By definition: we have m_D_E = m_D_E_marginal + m_D_E_infra;
m_D_E_marginal = sum(weight'.*(D_E_marginal_k));
m_D_E_infra    = sum(weight'.*(D_E_infra_k));
sd_theta       = (sum(weight'.*(theta_k-m_theta).^2))^0.5;
sd_D_E         = (sum(weight'.*((E_avg-E_0_avg)-m_D_E).^2))^0.5;
sd_D_E_marginal= (sum(weight'.*((D_E_marginal_k)-m_D_E_marginal).^2))^0.5;
sd_D_E_infra   = (sum(weight'.*((D_E_infra_k)-m_D_E_infra).^2))^0.5;
cov_theta_D_E  = sum(weight'.*((theta_k-m_theta).*((E_avg-E_0_avg)-m_D_E)));
corr_theta_D_E = cov_theta_D_E/(sd_theta*sd_D_E) ;
cov_theta_D_E_marginal  = sum(weight'.*((theta_k-m_theta).*((D_E_marginal_k)-m_D_E_marginal)));
corr_theta_D_E_marginal = cov_theta_D_E_marginal/(sd_theta*sd_D_E_marginal);
cov_theta_D_E_infra     = sum(weight'.*((theta_k-m_theta).*((D_E_infra_k)-m_D_E_infra)));
corr_theta_D_E_infra    = cov_theta_D_E_infra/(sd_theta*sd_D_E_infra);

%Reconstruct the welfare component based on sufficient statistics
% [yes, they are correct small difference of less than 1%]
D_Mis_SumStats_1 = mean(elec)*(m_theta-1)*m_D_E;
D_Mis_SumStats_2 = mean(elec)*corr_theta_D_E*sd_D_E*sd_theta;

%Compare whether the sufficient stats approach matches: D_Mis_SumStats vs D_Mis
D_Mis_SumStats = D_Mis_SumStats_1 + D_Mis_SumStats_2; 

%Compare the sufficient stats approach for the heterogeneity components
%for the infra vs marginal models.
D_Mis_SumStats_1_marginal = mean(elec)*(m_theta-1)*m_D_E_marginal;
D_Mis_SumStats_2_marginal = mean(elec)*corr_theta_D_E_marginal*sd_D_E_marginal*sd_theta;
D_Mis_SumStats_1_infra    = mean(elec)*(m_theta-1)*m_D_E_infra;
D_Mis_SumStats_2_infra    = mean(elec)*corr_theta_D_E_infra*sd_D_E_infra*sd_theta;

%Check that these calculations are correct: D_Mis_SumStats_check =
%D_Mis_SumStats [yes, they are correct]
D_Mis_SumStats_check = D_Mis_SumStats_1_marginal + D_Mis_SumStats_2_marginal + D_Mis_SumStats_1_infra + D_Mis_SumStats_2_infra;

end
